#include <assert.h>
#include "remez.h"
#include "remez_invsqrt.h"

using namespace mpfr;

// 要逼近的函数 f(x)= 1/x
mpreal invsqrt_func(mpreal x)
{
    if (x == mpreal(0) )
        return const_infinity(1); // 返回正无穷大
    
    return mpreal(1.0) / sqrt(x); // 返回1/x
}


// f(x)= 1.0/sqrt(x), 返回 f′(x)
mpreal invsqrt_derivative(mpreal x)
{
    return mpreal(-1.0) /( mpreal(2.0) * x * sqrt(x));  
}

// f(x)= 1/x
// my_invsqrt(x)  = coeffs[0] + coeffs[1]*x
// 返回 my_invsqrt(x) - f(x)的导数
// my_invsqrt′(x)- f′(x)
mpreal invsqrt_error_derivative(const std::vector<mpreal> &coeffs, mpreal x0,  int degree) // degree为多项式的阶数
{
    assert(degree+1<=coeffs.size());
    // coeffs=[c0,c1]
    // f(x)= c0 + c1*x
    // f′(x)=c1
    mpreal r= polynomial_derivative( coeffs, x0, degree);
    return r-invsqrt_derivative(x0);
}

// 计算 my_invsqrt(x)-invsqrt_func(x),即计算　my_invsqrt 在x点的误差
mpreal invsqrt_get_error_at_x( const std::vector<mpreal> &coeffs, mpreal x, int degree) // degree为多项式的阶数
{
    mpreal r = evaluatePolynomial(coeffs,x,degree);
    return r- invsqrt_func(x);
}
